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We present here two new algorithms for point-spread func- 
tion reconstruction that aim at suppressing the use of these Uij 
functions to avoid the storage of a large amount of data and to 
shorten the computation time of this PSF reconstruction. Both 
algorithms take advantage of the eigen decomposition of the 
residual parallel phase covariance matrix. In the first algorithm, 
the use of a basis in which the latter matrix is diagonal reduces 
the number of Uij functions to the number of mirror modes. 
In the second algorithm, this eigen decomposition is used to 
compute phase screens that follow the same statistics as the 
residual parallel phase covariance matrix, and thus suppress 
the need for these Uij functions. Our algorithms dramatically 
reduce the number of Uij functions to be computed for the 
point-spread function reconstruction. Adaptive optics simula- 
tions show the good accuracy of both algorithms to reconstruct 
the point-spread function. 
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Abstract. The knowledge of the point-spread function compensated by adaptive optics is of prime importance in several image 
restoration techniques such as deconvolution and astrometric/photometric algorithms. Wavefront-related data from the adaptive 
optics real-time computer can be used to accurately estimate the point-spread function in adaptive optics observations. The only 
point-spread function reconstruction algorithm implemented on astronomical adaptive optics system makes use of particular 
functions, named Uij. These Uij functions are derived from the mirror modes, and their number is proportional to the square 
number of these mirror modes. 



1. Introduction 

Since the advent of adaptive optics (AO), it has been possible 
to perform imaging with a spatial resolution very close to the 
diffraction Umit. Despite this large improvement, the correction 
of the atmospheric turbulence by an AO system is only partial, 
and point source images are for example still affected by a halo 
that surrounds them. 

Such effects of the partial AO correction can be corrected 
by image restoration techniques such as deconvolution algo- 
rithms. Except in the "blind" deconvolution case (Fusco et 
al. 1999; Christou et al. 1999), these algorithms need an es- 
timation of the point-spread function (PSF), or its correspond- 
ing optical transfer function (OTF) derived from the PSF by 
a Fourier transform, to restore the image unaffected by the at- 
mospheric turbulence. In addition, astrometric and photometric 
algorithms, e.g., StarFinder (Diolaiti et al. 2000) or DAOPHOT 
(Stetson 1987), usually also need an estimation of the PSF. 

Wavefront-related data delivered by the AO real-time com- 
puter can help to accurately estimate the PSF. Veran et al. 
(1997) have been the first to develop such a PSF reconstruc- 
tion algorithm. Implemented on the CFHT curvature sensing 
AO system PUEO (Arsenault et al. 1994), it has been rou- 
tinely delivering reconstructed on-axis PSFs for more than 8 
years now. A first attempt to transpose this algorithm to the 
Shack-Hartmaim (SH) wavefront sensor has been undertaken 
by Harder & Chelli (2000). 

Based on the Veran et al. (1997) algorithm, PSF reconstruc- 
tion has been developed for three AO systems, equipped this 
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time with SH wavefront sensors, and tested during a few runs 
of observations, leading to good results 

- Weiss (2003) has written a piece of PSF reconstruction soft- 
ware for ALFA (Kasper et al. 2000), the SH AO system of 
the Calar Alto 3.5m telescope; 

- Jolissaint et al. (2004) has written OPERA, a piece of PSF 
reconstruction software for Altair (Herriot et al. 2000), the 
4-quadrant SH AO system of the Gemini North telescope; 
and 

- Fitzgerald et al. (2004) have written a piece of PSF recon- 
struction software for the SH AO system of the UCO/Lick 
Observatorys 3 m Shane Telescope (Bauman et al. 2002). 

These current existing algorithms use particular functions, 
usually named Uij, computed from the mirror modes. The num- 
ber of Uij functions is proportional to the square number of 
mirror modes, which leads to gigabytes of data to handle for 
systems with 150 to 200 actuators, and thus limits the efficiency 
of the PSF reconstruction process. 

The goal of this article is to propose two new algorithms 
that avoid the use of these Uij functions. The global scheme of 
the PSF reconstruction is not affected; we have simply replaced 
the steps involving the Uij functions. We also provide, as a by- 
product, a way to estimate the PSF likelihood. The latter may 
be crucial for image deconvolution algorithms, which currently 
lack this kind of information. 

We present the main points of the PSF reconstruction al- 
gorithm developed by Veran et al. (1997), as well as their dif- 
ferent assumptions in Sect. 2. In Sect. 3, we describe the two 
algorithms we propose, and in Sect. 4, the tests from AO simu- 
lations. We discuss the choice and use of each algorithm in the 
last section. 
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2. The current Uij algorithm From the expression of the structure function of the residual 

« V , «^ . r^^^ parallel phase: 

2.1. The long-exposure AO-corrected PSF expression 

In the PSF reconstruction algorithm developed by Veran et D^^ (x,p) = (y>f^^(x) - <pq(x + p)j \ (4) 
al. (1997), assuming that the phase structure function defined 

in Eq. 4 below is stationary over the pupil, the AO-corrected and the decomposition of the phase on the basis of the mirror 

monochromatic long-exposure OTF is decomposed as follows: modes {M,(x)),=i n. 

10TF[pIA))=10TF^{pIA))x0TF,,{pIA), (1) A 

where: 



(5) 

!=1 



one obtains: 

N N 



- (oTF^^i^pl is the attenuation of the long-exposure OTF 

due to the partial correction of AO, and v v / \ 

- OT-F^P/a) is the perfect telescope OTF. ^S^^^P^ = ^ 2.<^ll'^ll^>KW-^'<^+^))KW-MXx+p)).(6) 



The phase 0^ can be split into two parts: which be- ^, ^ . r , ■ , , „ , , 

, ^ ^, ^ J u »u ■ ™ J J The mean structure function of the residual parallel phase 

longs to the vector space spanned by the mirror modes, and _ t- t- 

, u- u ■ ,u ^ , ,u f Da (p) is the mean of Z)<* (x, p) over jc: 

(pfj^, which IS orthogonal to the former space. '^'ii ^^w 

/ / x\ / / x\ / . x\ I D. (x,p)Pix)Pix-t-p)6x 

[otf{pIa)) = {otf,^(^pIa)) X {otf,^XpI^)) \(p) = ^ — I . (7) 



xOTF^{pIA) (2) J 



P{x) P{x -\- p)Ax 



and, this expression can be written The expression of D^^^(p) is similarly derived from 

D^^^{x,p). 

I , Y _ 1 _ The corresponding AO-corrected monochromatic long ex- 

10TF{p/a)\ * exp ( - -D^^^ (p)) X exp ( - -D^^^ (p)) posure PSF is derived as the Fourier transform of the OTF 



X OTFta{p/A), (3) 



where: 



- (oTF^^^ ipl'^ attenuation of the long-exposure OTF 



2.2. The Uijifi) functions 
Equation 6 can be rewritten: 

N N 



The Uijip) functions are defined by: 



due the mirror component of the phase, i.e., the "residual _ V V .^h^.^ Uijip) (8) 

parallel phase", ^" 

- (oT F ^^ {p I A^ is the attenuation of the long-exposure OTF 

due the component of the phase belonging to the space 

perpendicular to the mirror space, i.e., the "perpendicular r 

phase", J (*^,- W - Miix + p)){Mj{x) - Mj{x + p))P{x)P{x + p)dx 

- D^^ (p) is the mean structure function of the residual paral- ^ >(9) 

lei phase, I P(x)P(x + p)6x 

- D^^ (p) is the mean structure function of the perpendicular 

phase, where Pir) is the pupil function and x a coordinate vector in 

- p is a pupil plane coordinate vector, and the pupil plane. 

- /I is the wavelength of observation. In practice, using the Fourier transform and the properties 



This decomposition is based on several assumptions (cf. 
Veran et al. 1997): 



of the correlation function (Veran et al. 1997), the Uij(p) func- 
tions are computed as: 



1 . the complex field amplitude is uniform over the pupil (scin- M 2!R (T{MiMjP)T* (P)- T{MiP)T{MjP)^ j 
filiation is neglected), f/r(p) = ^ L (lO) 

2. the residual parallel phase <;*q|(x) and the perpendicular T'^ {\ T{P) 1^) 
phase (peS^) follow Gaussian statistics, 

3. the phase structure function is stationary over the telescope where T is the Fourier transform function and % the complex 
pupil, and number real part function. 

4. the correlation between the mirror component and the per- Equation 8 is a key one for the experimental reconstruc- 
pendicular component of the phase is negligible. tion of PSFs. The covariance matrix (en en') has to be measured 
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experimentally on the AO system itself, by averaging the cross- 
products of wavefront measurements obtained during the time 
of the image exposure. 

In the current PSF reconstruction algorithms, derived from 
Veran et al. (1997), the matrix <e||e||') is the basic entry point 
from which one can deduce successively the phase structure 
function, the OTF, and then the PSF. Additionally, one has to 
compute, store once for all, and also read the Uijip) functions 
during the reconstruction process . 

In Eq. 8, the i and j indices play a symmetric role, so that 
there are actually Nx {N+ 1)/2 "useful" Uij{p) functions. As an 
example, in the case of the VLT AO system NAOS (Rousset et 
al. 2000), the 159 compensated modes lead to 12720 "useful" 
Uijip) functions. Today, the large number of Uij{p) hence rep- 
resents, depending on the array size and data type, several gi- 
gabytes of data to compute, store, and read. Even if in practice, 
Eq. 8 can be efficiently implemented by Fourier transform (cf. 
Eq. 10), this leads to a heavy PSF reconstruction process, which 
will turn out to be impossible to handle in the future since the 
next AO systems are expected to have a largely increased num- 
ber of modes: about 1370 actuators for the VLT-Planet Finder 
AO system (Fusco et al. 2005) and several tens of thousands for 
extremely large telescopes. In the following, we propose a way 
to achieve this computation, starting from the same covariance 
matrix {e\\e\\), without using the Uijip). 

3. Theory of the proposed algorithms 

3.1. The Vii algorithm 

Let us consider the vector e\\{t), hereafter €\\, made of the 
{^111)1=1...^ coefficients, i.e., the vector representing 0q|(ji:, in 
the basis of the mirror modes M, (ac). The eigen decomposition 
of the residual parallel phase covariance matrix is: 



A = B'{€«e{)B, 



(11) 



where A is a diagonal matrix that contains the [/ii\i=\,„N eigen- 
values and B is the matrix of eigenvectors: B'B = BB' = Id. 
Equation 1 1 can be written: 



A = ((B'e||)(B'e||)'). 



(12) 



The vector rj equal to B'€\\ represents (l)^^{x,t) in ffie basis 
that diagonalizes the residual parallel phase covariance matrix. 
Its coefficients are noted {)7,},=i...Ar. From Eq. 12, the covariance 
matrix {qif) is diagonal, i.e., in this new basis, the residual par- 
allel phase covariance matrix is diagonal, and the mean residual 
parallel phase structure function reduces to: 



(13) 



i=\ 



where the Vijip) functions are ffie equivalent in the new basis 
to the Uijip) functions (Eq. 8). Similarly to Eq. 9, the Vijip) 
functions are defined by 



^ [M\ix) - M[ix + p))[M'.ix) - M'jix + p))Pix)Pix + p)dx 



such ffiat M', the matrix made of the eigenvector modes 
{M'.ix)]i=\,„N is given by M' = B'M, M being ffie matrix made 
of the mirror modes {Miix)}i=\ ^. 

Using ffie V,, algorithm, ffie computation of ffie residual 
parallel phase OTF only requires the computation of a number 
A' of functions Vuip). However, these Vuip) functions have to 
be computed on the fly for each estimation of ffie mean residual 
parallel phase structure function. 

3.2. The "instantaneous PSF" algorithm 

The solution we propose here is similar to the algorithm pre- 
sented by Roddier (1990) to simulate atmospherically distorted 
wavefronts, calculated from the covariance matrix of the coef- 
ficients of their expansion in Zernike modes. We extend this 
algorithm to any modal basis and covariance matrix, and use it 
to reconstruct AO-corrected PSFs. 

Let us consider again the eigen decomposition of ffie resid- 
ual parallel phase covariance matrix: 



(eiieil') = BAB', 



(15) 



If one generates a vector rj whose coefficients are indepen- 
dent Gaussian random variables wiffi zero mean and variance 
equal to the eigenvalue A,, i.e., (rp]'} = A, ffien ffie vector 
j3 = Br] is a set of correlated random variables whose covari- 
ance matrix is {e\\e\\'y. 



(jSp'} = {Br^iB') = BAB' = {e^^q{). 

The phase represented by the vector p is: 



(16) 



,f>ix,t) = YjMt)Miix), 



(17) 



1=1 



and ffie "instantaneous" PSF corresponding to that phase is: 



PSF\\ix,t) = ||:r(exp(j0(x,O 



(18) 



Then, by generating random rj vectors such that {rjri'} = A, 
we build instantaneous PSFs ffiat, on average, converge to ffie 
long-exposure PSF of ffie mirror space. Note ffiat the latter is 
not the "full PSF" that would be observed at the telescope since 
it does not include the uncorrected part of ffie phase (cf. Eq. 2). 
Finally: 

(OTF^^^ (x/^)) X OTFt,i{x/A) =rn] PSF^^ix, t)\ . (19) 



4. Test of the algorithms 

4.1. Description of the simulation 



I 



Fix) Fix + p)6x 



To test our algorithms, we have used a Monte Carlo-based AO 
simulation software developed at ONERA. It is a complete end- 
to-end AO simulator divided into four main modules (calibra- 
tion, propagation, closed-loop, and focal plane imaging mod- 
''-^^-'ules) ffiat make it as close as possible to an actual system. The 
complete algorithm is fully described in Conan R. et al. (2004) 
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Fig. 1. Strehl ratio at 2.2 fim vs. guide star visual magnitude for a 
NAOS-like AO system. 

and has been used extensively for the design of AO systems, 
especially the future Planet Finder System for the VLT (Fusco 
et al. 2005), as well as for the tests and vahdations of existing 
AO systems (NAOS for instance). 

We have tested our algorithms in the simple case of a 
NAOS-like AO system (Rousset et al. 2000): a 14x14 sub- 
pupil Shack-Hartmann wavefront sensor with 8x8 pixels per 
sub-aperture, a read-out noise of 3 e" per pixel, and a sampling 
frequency of 500 Hz, 2048 loop iterations. The correction was 
performed with a 185 actuator piezostack deformable mirror 
plus a tip-tilt mirror. The wavefront sensing and observation 
wavelengths of the simulation were 0.65 and 2.2 jum, respec- 
tively. The seeing was 0.85" at 0.5 fim. 

Since we aimed at testing our algorithms with different con- 
ditions of correction, we ran the simulation with a guide star 
magnitude ranging from 7 to 15 so that the resulting Strehl ra- 
tio ranged from ~7Q% down to ~0.2% (Fig. 1). For a given 
guide star magnitude, we stored all the values e\\(t) obtained 
from the simulation to compute the covariance matrix (eyeii') 
and ran the Uij, V,,, and "instantaneous PSF" algorithms (code 
detailed in Appendix A) from this covariance matrix to derive 
the corresponding OTPs for comparison. 

4.2. Results of the simulations 

Figure 2 represents the circular mean of the residual "atmo- 
spheric" OTFs (i.e., not multipUed by the telescope OTF) pro- 
duced by the different algorithms in three cases: good, mod- 
erate and poor correction (7.3"', 12.3"', and 13.6"" magnitude 
guide stars, respectively). Within the numerical uncertainties, 
the Uij and V,, algorithms produce exactly the same OTFs, 
which is not unexpected since they mathematically do the same 
computations, but in a different modal basis. Their correspond- 
ing PSF profiles in Fig. 3 are consequently also the same. 

This is of course not the case for the "instantaneous PSF" 
algorithm, which needs a large enough number of iterations to 
converge. However, even in the poor correction case, the OTF 
profile obtained with the "instantaneous PSF" algorithm well 
reproduces the Uij profile (Fig. 2) with a maximum error of 
a few 10"^ around 0.05 D/A. This is also confirmed with the 




Normalfzed frequency (D/X units) 

Fig. 2. Circular mean of the "atmospheric" OTFs for different condi- 
tions of correction: from top to bottom, the guide star magnitude is 
7.3, 12.3, and 13.6. For each guide star magnitude, the solid line is 
obtained with the Uij/Vu and the crosses with the "instantaneous PSF" 
algorithm". 
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Fig. 3. Circular mean of the PSFs for different conditions of correc- 
tion: from top to bottom, the guide star magnitude is 13.6, 12.3, and 
7.3. For each guide star magnitude, the solid line is obtained with the 
Uij/Vii and the crosses with the "instantaneous PSF" algorithm". Note 
that the PSFs considered here only correspond to the residuals and do 
not include the perpendicular part of the phase. 

PSF profiles that are almost superimposed, even in the poor 
correction case (Fig. 3). 

5. Discussion 

The Uij and Vu algorithms mathematically produce the same 
OTFs. In practice, the Uij algorithm requires reading each 
of the A' X (A^ -I- l)/2 Uij functions (where N is the number 
of modes), which have been computed and stored previously, 
once. In comparison, the V,, algorithm requires to diagonahze 
the covariance matrix, to compute the new modes [M'.(x)]i=\,„N, 
and to compute each Vu function. This latter step is composed 
of calls to these functions (cf. Appendix A): 

- 2FFTs, 

- 1 square function, 

- 1 real part of a complex number, and 
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- 1 modulus of a complex number. 

As an example, the V,, algorithm took ~8 seconds in our 
simulation: ~2 seconds to diagonaUze the covariance matrix 
and compute the modes (the former taking a negligible time), 
and ~6 seconds to compute the A' V,, functions. This is ~25 
times faster than the Uij algorithm (~191 seconds). In addition 
to this huge gain in computation time, a large amount of disk 
space is saved. This causes the Vu algorithm to always be pre- 
ferred to the Uij one. 

As noticed by Conan J.-M. (1994), averaging short- 
exposure OTPs, as we do in the "instantaneous PSF" algorithm, 
is a process that converges very slowly, especially at large D/tq 
or a low correction level. In addition, it does not lead to the in- 
finitely long exposure OTF since a given number of short ex- 
posures are averaged. Besides, Conan J.-M. (1994) has shown 
that in the poor correction case, the error in computing the long- 
exposure OTF of such an algorithm is larger than for the Uij 
algorithm, and then the V,, algorithm as well. 

Though, we emphasize that, in addition to the OTF compu- 
tation itself, the "instantaneous PSF" algorithm can provide an 
estimation about the variabiUty of the OTF, which can help a 
lot in some deconvolution algorithms. The estimation of the in- 
finitely long exposure OTF that results from the convergence 
of our "instantaneous PSF" algorithm and that corresponds 
to a given covariance matrix is unique: let us call it 

OTFooip/A). The dispersion in the random generation of OTPs 
can be computed as 

a^(p/A) - (^\OTF^(p/A) - OTFi(p/A)\f^^, (20) 

where OTFi is the draw of a randomly-generated OTF. If we 
call OTFohs(p/A) the OTF actually observed on the instrument 
during the given, non infinite, observing time Tint, when the 
given covariance matrix was measured, we can evaluate 

how far our estimation OTF oo is from OTFobs by writing: 

(]\OTF^(p/A) - OTFobs(p/4f) = (21) 

where n is the equivalent number of independent realisations 
of PSFs, whose sum has resulted in the final PSF observed by 
the instrument during the given, non infinite integration time 
Tint. An estimation of n could be obtained for example from a 
full simulation of the AO system under the same atmospheric 
conditions as during the observation. It is also be reasonable to 
consider that the impact of the correction by the AO system is 
to shorten the image correlation time compared to the image 
correlation time to(A) of the atmospheric seeing (Rigaut et al. 
1991). We can then find a lower bound given by « > Tiat/To(A). 
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Appendix A: IDL codes 

A.1. Common parts 

cbmes = readfits('cbmes.fits') 
s2m = readfits('slopes2modes. fits') 
modes = readfits('modalbasis.fits') 
nmodes = (size(modes)) (3 ) 
lambda jm=2.2d-6 
lambda^so=0.65d-6 
ratio Jambda=lambda^so/lambda_im 
; Telescope OTF computation and corresponding mask 
apert=readfits('pupille.fits') 
dim= (size(apert))( 1 ) 
dim2=2*dim 
pup=dblarr(dim2,dim2) 
pup(0,0)=apert 
pupfft=fi^(double(pup)) 
conjpupflit=conj(pupfft) 

otftel = real_part(fft(pupfft*conjpupff't,/inverse)) 
den=l/(otftel) 

den(where(finite(den) eq 0))=0 
mask=fltarr(dim2,dim2)-i- 1 
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mask(where(otftellt le-5))=0 > 
otftel=otftel/max(otftel) 
;Computation of the covariance matrix over the modes 
cbmode = cbmes##transpose(s2m) 
nbcb - (size(cbmode))(2) 
covmode = (cbmode#transpose(cbmode))/nbcb 

A.2. The Uij algorithm 

; Computation and storage of the Uij functions 
modei = dblarr(dim2,dim2) 
modej = dblarr(dim2,dim2) 
for i=0,nmodes-l do begin 

modei(0,0) = modes(*,*,i)*apert 
for j=i,nmodes-l do begin 

modej(0,0) = modes(*,*,j)*apert 
terml - real_part(fft(modei*modej)*conjpupfft) 
term2 = real_part(fft(modei)*conj(fft(modej))) 
uij = real_part(fFt(2*(terml-term2),/inverse))*den*mask 
writefits,'Uij/u_' +string(i,format=' (i2.2)')+' -' +$ 
string(j,format=' (i2.2)')+' .fits' ,uij 

endfor 
endfor 

; Computation of the OTF with the Uij 
dphl - dblarr(dim2,dim2) 
for i=0,nmodes-l do begin 
for j=i,nmodes-l do begin 

uij =readfits('Uij/u_'+string(i,format='(i2.2)')+'_'+$ 

string(j ,format= ' (i2.2) ' )+ ' .fits ' ,/silent) 
fac - double((ine j)+l) 
dphl = dphl+(fac*covmode(i,j)*uij) 
endfor 
endfor 

otfl - exp(-0.5*dphl*ratioJambda''2)*mask 
otfl = otfl/max(otfl) 

A.3. The Vii algorithm 

; New modes that diagonaUze the covariance matrix 
1 = (eigenql(covmode,eigenvectors=b))>0 
s = b#diag_matrix(sqrt(l)) 

newmodes = reform((reform(modes,dim''2,rmiodes))#b,$ 

dim,dim,rmiodes) 

; Computation of the OTF with the Vii 
newmodei = dblarr(dim2,dim2) 
temp - dblarr(dim2,dim2) 
for i=0,nmodes-l do begin 

newmodei(0,0) = newmodes(*,*,i)*apert 

terml = reaLpart(fFt(newmodei''2)*conjpupfft) 

term2 = (abs(fft(newmodei)))"2 

temp = temp+((terml-term2)*l(i)) 
endfor 

otf2 = exp(-0.5*real4)art(fFt(2*temp,/inverse))*den$ 

*ratioJambda"2)*mask 

otf2 = otf2/max(otf2) 



.4. The "instantaneous PSF" algorithm 
psf = dblarr(dim2,dim2) 

tmp = ratioJambda*reform(modes,dim''2,nmodes) 

for i=0,nbcb-l do begin 

phi(0,0) = reform((s#randomn(seed,nmodes))##tmp,dim,dim) 
psf = psf+(abs(fft(pup*exp(dcomplex(0,phi)))))'"2 

endfor 

psf = psf/nbcb 

otO = real_part(fft(psf)) 

otf3 = otf3/max(otf3)/otftel*mask 
otf3(where(finite(otf3) eq 0))=0 



